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to turducken (turduckens, turduckening, turduckened, turduckened) [math.]: To stuff a black hole. 



We analyze and apply an alternative to black hole excision based on smoothing the interior of 
black holes with arbitrary - possibly constraint violating - initial data, and solving the vacuum 
Einstein evolution equations everywhere. By deriving the constraint propagation system for our 
hyperbolic formulation of the BSSN evolution system we rigorously prove that the constraints prop- 
agate causally and so any constraint violations introduced inside the black holes cannot affect the 
exterior spacetime. (This does not follow from the causal structure of the spacetime as is often 
assumed.) We present numerical evolutions of Cook-Pfeiffer binary black hole initial configurations 
showing that these techniques appear to work robustly for generic data. We also present numerical 
evidence from spherically symmetric evolutions that for the gauge conditions used the same station- 
ary end-state is approached irrespective of the choice of initial data and smoothing procedure. 

PACS numbers: 04.20.-q,04.25.Dm,04.30.Db 



I. INTRODUCTION 

Currently, there are essentially two different ways of 
dealing with singularities in the numerical evolution of 
orbiting black holes. One technique is black hole excision 
[H, 0], where the interior of each black hole is removed 
from the computational domain by an inner boundary. 
The other is the so called "moving punctures" technique 
@j El j where the initial asymptotically flat regions inside 
each black hole are represented by "puncture points" . 
Long, multi-orbit binary black hole simulations have been 
achieved over the last few years using both excision 0, Q 
and moving punctures (see Q and references therein). 

The puncture technique does not make use of the black 
hole excision idea, at least not in the classical sense 
of placing an inner boundary inside each black hole. 
Instead, the fields that initially describe the puncture 
points are allowed to evolve freely in the (topologically 
trivial) computational domain and the subtleties of black 
hole excision are replaced by the subtleties involved in 
approximating the singularities in the equations at the 
puncture points. The particular appeal of the "moving 
punctures" technique compared to black hole excision is 
that it appears to be simpler to achieve a stable dis- 
cretization near the puncture points than at an excision 
boundary; however, there appears to be an implicit limi- 
tation of the method in that it is in principle tied to the 
use of puncture data. Recently, light has been shed on 
the geometric picture behind moving punctures [1, H, |T(| • 

In this paper we discuss a technique for evolving black 
holes which shares the simplicity of moving punctures 
but is not restricted to puncture-type initial data and 



does not need any regularization of the equations near 
special points. The method also relies on the intuitive 
idea behind black hole excision that "no physical infor- 
mation can escape from the interior of a black hole" , but 
proceeds in a different way. In particular, it does not 
require placing an inner boundary per black hole in or- 
der to remove the interiors. The computational domain 
in this technique is trivial (from a topological point of 
view) and the discretization therefore remains simple. 

The basic idea is the following: if no physical infor- 
mation can leave the interior of the black hole, why not 
just change the interior to one's advantage? The spirit of 
this idea is not new, and has been advocated for a long 
time in several forms, most notably by Bona and collab- 
orators [ill, Ell EH and by Misner [3]. In particular, in 
fl2| . a "free black evolution" approach was advocated, 
where the interior of each black hole is smoothed with 
arbitrary data and the vacuum Einstein evolution equa- 
tions are solved everywhere. In general the smoothing 
process generates constraint violations. Thus, a key in- 
gredient of this approach is to guarantee that the form of 
the equations does not allow for constraint violations to 
propagate to the outside. This is highly non-trivial. In 
fact, it is well known that depending on the form of the 
Einstein equations used, gauge and constraint modes can 
propagate with arbitrary (including supcrluminal) speeds 
and, in particular, constraint violations can leak from the 
interior of black holes to the outside. Though we use a 
different formulation of the equations (a version of BSSN 
as opposed to the Z4 system [H|]), this "free black hole 
evolution" approach is exactly the one that we analyze 
and apply in this paper. Even though in several aspects 
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this is different from the "stuffed black hole proposal" 
[ill ], we will refer to our particular implementation as 
the relativistic turducken Ha . 



II. NO CONSTRAINT LEAKING IN THE 
TURDUCKENING: AN ANALYTICAL PROOF 

The exact version of the BSSN system we use is given 
by Eqs. (3)-(7), (21)-(23) in [l6| where we set the param- 
eter m to one and the source terms S, Sij and S l to zero. 
Furthermore, the lapse a and the shift f3 l are evolved 
according to the 1 + log slicing condition boa = —2aK 
and the "hyperbolic Gamma driver" [l7| like conditions 
<9 /3 4 = 3BV4, d B l = dot 1 - B l /2 , respectively, where 
do = dt — ftdj. The term doY % in the last equation is set 
equal to the right-hand side of the evolution equations for 
the f 1 symbols. As noted in [lj[ the use of do (as opposed 
to dt) in the above equations simplifies the analysis of the 
hyperbolic structure of the equations. Later, it was also 
found to be important in practice for long-term binary 
evolutions In addition, using do for the lapse implies 
that the slicing obtained is independent of the choice of 
shift vector @. 

The wcll-posedness of the resulting Cauchy problem 
was analyzed in [l(|. A sufficient condition for well- 
posedness is strong hyperbolicity of the evolution equa- 
tions. (See [TjJ Ho[ for definitions that apply to second 
order systems.) In our case, the equations are strongly 
hyperbolic if and only if the lapse a and the conformal 
factor 4> are smooth functions satisfying a > 0, \<f>\ < 00 
and h := 2a — e 4 ^ 7^ 0. The last condition is typically 
violated, at least on some two-surface. This is so because 
in general, a — -> 1 and — > and therefore h — > 1 as one 
approaches the main asymptotically flat end, while near 
black holes a is small and <f> is large and positive (for the 
coordinate conditions used here typically a « 0.3 at the 
horizon, and a — > and ^ 00 at any punctures) so that 
h < near a horizon. Therefore, the function h must be 
zero somewhere in between. On the other hand, if the re- 
gions where h = are, for example, sets of zero-measure 
in the computational domain there is hope that the vio- 
lation of the condition h ^ still allows for a well posed 
Cauchy problem. The numerical simulations in Sect. IIVI 
below show no apparent sign of numerical instability. 

The characteristic speeds (with respect to normal ob- 
servers) for our evolution equations are the following 
0,±l,±/ii,±^ 2 ,±M3, where /ii = y/2/a, ^2 = 
V3e 2(t> /2a, ^ = e^/a. It is possible to give a precise 
meaning to the different characteristic fields and speeds 
in the high-frequency limit [2ll, |22|. In that limit, fields 
propagating with speeds fi\, /J2 and /13 correspond to 
gauge modes, while the fields corresponding to gravi- 
tational radiation and constraint-violating modes have 
speeds ±1 and 0, ±1 respectively. As we will see below, 
the constraint propagation system possesses the charac- 
teristic speeds and ±1. 



The BSSN system is subject to the Hamiltonian and 
momentum constraints H = and Mi = plus three 
extra constraints associated with the introduction of 
the F l symbol as independent variables, namely Cp := 
r 4 + dj^ = 0, where 7^ refers to the inverse of the con- 
formal metric. In order to obtain a solution to Einstein's 
vacuum field equations, these constraints have to be sat- 
isfied. We now show that it is sufficient to solve them 
on an initial Cauchy surface in the region exterior to 
the black holes. The constraint propagation system then 
guarantees 1 that these constraints hold at every time fu- 
ture to the initial surface and at every point outside the 
black hole regions, independent of any constraint viola- 
tion in the interior of the black holes. We show this by 
deriving the constraint propagation system and casting 
it into first order symmetric hyperbolic form. Then the 
causal propagation of the constraints can be shown via a 
standard energy inequality provided all the characteristic 
speeds (as measured by normal observers) of the system 
are smaller than or equal to one in magnitude. 

Using the Bianchi identities, imposing the evolution 
equations and introducing the additional constraint vari- 
ables Zij = (diCp)'jkj, the constraint propagation system 
can be rewritten as a first order system of the form 

doC = a [A( u yd l C + B(u)C] , (1) 

where C are the constraint variables, u — 
(a, 4>, "fij, K, Aij, r l ) are the main variables, and 
A 1 , A 2 , A 3 and B are matrix- valued functions of u. De- 
composing Zij = Z^ij) + Zuj] + jijZ/S into its trace-free 
symmetric part Zuj)-, its anti-symmetric part Zuj], and 
its trace Z = j lJ Z^ with respect to the physical three- 
metric 7,j , and representing C in terms of the variables 
C = (C r , Si := 2mH+Z, S 2 := H+2aZ, M j} %),%]), 
the principal symbol A(n) = A.(u) l ni is given by 

A(n)C = (0, 0, n j Mj , ± nj S 2 + \n l Z {l]) + \n'Z M , 

2(n {l M j) ) TF ,2n [t M ]] y (2) 

where ri 1 = ~l lJ nj and is normalized such that n^n 4 = 
1. This system is symmetric hyperbolic, and its charac- 
teristic speeds (with respect to normal observers) are 
and ±1. A symmetrizer is given by the quadratic form 

C T HC = Tij-CfCf + S 2 1 +\sl+"f i] M l M j 

+ \i lk i ]l z^)Z (kl) + \i lk ^ l z m z [kl] . 

The symmetrizer, along with the fact that there are no 
super luminal characteristic speeds, allow us to obtain an 



1 However, constraint violations can still be introduced by im- 
proper outer boundary conditions. 
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energy estimate for the constraint variables C and to 
show that no constraint violations from the interior of 
a black hole can propagate to the outside. The explicit 
estimate will be presented elsewhere along with more de- 
tails of the results presented in this paper. 



III. SINGLE BLACK HOLE EVOLUTIONS AND 
THE END STATE 

In this section we present insights obtained by applying 
the turducken technique to a single spherically symmetric 
black hole. For these studies we use both the three- 
dimensional (3D) code described in section IV, as well as 
the one dimensional (ID) BSSN code discussed in [23| . 
Both codes use a formulation of the BSSN equations that 
is strongly hyperbolic everywhere except in regions of 
the computational domain that are likely sets of measure 
zero, and have causal constraint propagation. 

For a single black hole we use turduckened Kcrr-Schild 
(KS) initial data. Without turduckening, a KS slice hits 
the singularity. We first define the spacetime metric 
flV = + ^H£^£ u in terms of Cartesian coordinates 
x, y, z, where r]^ is the Minkowski metric, H = M/f 
and ifj, = (1, x, y, z)/f. Here, f is defined in terms of co- 
ordinate radius r = (x 2 + y 2 + z 2 ) 1 ! 2 by f = (r p + e p ) 1 / p . 
The contravariant metric g^ u is obtained from g^ by 
raising indices with rj^ . In Cartesian coordinates the 
initial metric is defined by the spatial components of g^ 
and the initial extrinsic curvature is defined by the usual 
expression K vi = (-g lj +(3 k d k g i: j+2g k ( l d j )l3 k )/2a, where 

a = 1/ \J-g tt and ft = -g lt /g a . The initial data for the 
ID code is obtained by transforming the Cartesian data 
to spherical coordinates. 

For r » t we find r sa r and the initial data coincides 
with a KS slice of a non-rotating black hole. For r close 
to the origin, the data are smooth and regular as long 
as e 0. This form of turduckening is not ideal since 
it leads to constraint violations that extend beyond the 
horizon r — 2 M. Typical values used in our simulations 
are e = 0.1M and p = 4. These values lead to initial 
violations of the Hamiltonian constraint of ~ 10 4 /M 2 at 
r = and - 10~ 5 /M 2 at r = 2 M. 

Experiments in ID show that after an evolution time 
of 50 M, the Hamiltonian constraint violation throughout 
the computational domain drops to a level ~ 10 _5 /M 2 . 
Similar results hold for the other constraints. By t = 
50 M the data have become nearly stationary; the final 
state in the t — > oo limit coincides with a portion of the 
stationary 1 + log slice of Schwarzschild. This is the same 
end state obtained with puncture evolution [1, [l(| ■ The 
key ingredient responsible for these remarkable behaviors 
is the Gamma-driver shift condition. With this condition 
the shift grows large in the interior region to counteract 
the grid stretching that would otherwise occur as the 
lapse collapses. As a result the time flow vector field tips 
outside the physical light cone (toward increasing r) and 
the grid points near r = are quickly driven out of causal 
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Figure 1: Areal radius R versus proper distance d from the 
horizon. The initially turduckened KS become indistinguish- 
able from the portion d >, — 9 of a stationary 1 + log slice after 
t « 50 M. The region R < 2M is the black hole interior. 

contact with the constraint violating portion of the initial 
data. With the constraints (nearly) satisfied everywhere 
in the computational domain, the numerical data repre- 
sents a slice of Schwarzschild that extends from region I 
of the Kruskal diagram, crosses the black hole horizon, 
and terminates at a resolution-dependent location inside 
the black hole. The 1 + log slicing condition then guides 
the slice to a stationary state. 

Fig. [T] shows the areal radius R versus proper distance 
d (in the radial direction) for a single non-spinning black 
hole, obtained from the ID code with resolution M/200. 
The data evolve to the stationary 1 + log slice in spite of 
the fact that the initial data violate the constraints. 

Alternatively, the initial KS data can be changed only 
inside a sphere of radius r < 2A/. In ID simulations such 
initial data can lead to the formation of gauge shocks, like 
those discussed in Q . The shocks typically form just out- 
side the black hole, independent of the parameter ro; this 
suggests that the formation of shocks is a consequence of 
the gauge conditions, and not the black hole turducken- 
ing. We have not seen this behavior in 3D, perhaps due 
to lack of resolution. 



IV. BINARY BLACK HOLE EVOLUTIONS 
USING COOK-PFEIFFER DATA 

We evolve quasi-equilibrium binary black hole initial 
data using the form of the equations described above, im- 
plemented in CCATIE, a three-dimensional adaptive mesh 
refinement code which uses the Cactus framework [24[ 
and the Carpet mesh refinement driver [25|, [2(|. This 
evolution code is fourth order accurate. It uses cen- 
tered finite differencing operators, except for the advec- 
tion terms which are upwinded. We use fifth order spa- 
tial and second order temporal polynomial interpolation 
at mesh refinement boundaries, and buffer zones as de- 



4 



Black hole trajectories and shapes 
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Figure 2: Apparent horizons in a binary black hole evolution 
of Cook-Pfeiffer data. This figure shows the tracks of the 
centroids of the apparent horizons, as well as their shapes at 
t = M, at merger (t « 37 M) , and at late times (t > 200 M). 



scribed in [25( to ensure stability. We therefore expect 
our code to be third order accurate in the limit of infinite 
resolution, and expect it to show approximate fourth or- 
der convergence away from the outer boundary and for 
the resolutions used here. We use a fourth order Rungc- 
Kutta time integrator with a CFL factor of 0.4. We use 
Sommerfeld outer boundary conditions for the individual 
components of the evolved variables, which are not con- 
straint preserving; we therefore place the outer bound- 
aries at a large distance from the source. 



The initial data were provided by Harald Pfeiffer [27 1 
and are described in [__| __|]. In particular we use 
the data set sep_07.00_59a.tgz in which the binary 
black hole system is expected to orbit approximately 
once before merging. These data have an ADM mass 
Madm ~ 2.44449 « 0.977795 M, where we use a scale 
factor M = 2.5. The black holes are centered about 
x = ±1.4 M, and the apparent horizons have a coordi- 
nate radius tah ~ 0.35 M. 

Our simulations use reflection symmetry about z = 
and 7r-rotation symmetry about the z axis. We choose 
a simulation domain with outer boundaries at 204.8 M. 
and use altogether 9 successively smaller levels of mesh 
refinement, where the finest level has an extent of 0.8 M, 
centered about each black hole. Our resolution is h = 
3.2 M on the coarsest grid, h = 0.8 M near the gravita- 
tional wave detector, and h = 0.0125 M on the finest grid. 
We include results from two coarser runs with coarse grid 
resolutions h sa 4.5 M and h w 4.1 M, respectively. 

The initial data are provided in terms of spectral ex- 
pansion coefficients for the ADM variables on multiple 
domains and need to be interpolated to our grid points. 
The initial data setup excises the apparent horizons but 
extrapolates a distance of up to 0.25 tah = 0.0875 M 
into the horizon. The remainder of the interior of the 
apparent horizons needs to be turduckened. 

We have experimented with various methods for tur- 
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Figure 3: Real part of the waveform Q^i extracted at -R = 
50 M. The vertical lines indicate approximately when the 
initial burst of spurious radiation first reaches the detector 
and when the common horizon is "seen" by the detector. The 
"junk" radiation near t — 50 M is a well-known feature from 
puncture evolutions. 



duckening the black hole interior, and find that the de- 
tails do not matter much in practice, as long as the space- 
time remains unchanged within the finite differencing 
stencil radius of the horizon. Since there are preciously 
few grid points between the excised region and the hori- 
zon, we chose a method which leaves all given spacetime 
data unchanged and fills in the excised points in a smooth 
manner. (One alternative would be a blending method 
which fills the excised region with arbitrary data, and 
then modifies some of the non-excised grid points to cre- 
ate a smooth match.) 

In particular, we solve the elliptic equation (d 6 /dx 6 + 
d 6 /dy 6 + d 6 /dz 6 )A = to fill the excised points of a 
quantity A, using standard centered derivatives every- 
where and using the given non-excised data as boundary 
conditions where necessary. This is equivalent to provid- 
ing boundary conditions for A and its normal derivatives 
dA/dn, and d 2 A/dn 2 . The result is therefore C 2 every- 
where within the horizon. We solve this equation with a 
standard conjugate gradient method. 

We follow the evolution of these data through merger 
and ringdown for about 200 M. Fig. __ shows the loca- 
tions, shapes, and tracks of the individual and the com- 
mon apparent horizons. A common horizon appears at 
about t = 37 M. The common horizon initially has a 
strong Y22 deformation which is radiated away. This is 
clearly shown in the real part of the I — m = 2 mode of 
the Zerilli function Q + , extracted on a coordinate sphere 
at R = 50 M and shown in Fig. __ Fig. __ shows the re- 
sults of a convergence test, although the resolutions are 
too close together to give reliable results. Both the hori- 
zon dynamics and waveforms are very similar to those 
from puncture initial data. We will present a study of 
this and other systems with larger initial separations in 
more detail elsewhere. 
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Waveform convergence 
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Figure 4: Difference between waveforms, scaled for 4th order 
convergence. The waveform phases have been shifted in time 
so that all resolutions have the same phase at the beginning 
of the merger radiation burst at t m 60 M. The initial high- 
frequency oscillations in the error are caused by the small 
amount of "junk" radiation near t = 50 M, and similar oscil- 
lations near t = 140 M are probably caused by its reflection 
at the coarsest refinement level. The noise near t — 170 M 
appears at a time when the waveform has already rung down. 



V. FINAL REMARKS 

A key property needed in a "free black hole evolu- 
tion" approach is that the constraints propagate causally. 
This cannot be taken for granted, and must be proved 
(or tested) for any particular formulation of the Einstein 
equations used. Note that even apparently small varia- 
tions in the evolution system can change the constraint 
propagation from causal to acausal. 

Causal propagation of the constraints alone is not suf- 
ficient. In modifying the initial data by smoothing away 
the singularity, we are not guaranteeing that the evolu- 
tion will proceed to a smooth, regular end-state. That 
this end-state is numerically well-behaved is the other 



key ingredient in any evolution that relics on modifying 
the interior of the horizon in some way. As the numer- 
ical evidence presented here shows, evolutions in spheri- 
cal symmetry do tend to a recognizable end-state for the 
given set of gauge conditions and form of the equations. 
It seems likely that a similar picture will hold away from 
spherical symmetry. 

Our work suggests that the turducken technique will 
hold irrespective of how and when the data inside the 
horizon are modified, thus allowing the method to be 
applied without modification to the final stages of evo- 
lutions performed with possibly different codes and/or 
methods, or to horizons formed e.g. in stellar collapse 
scenarios. 

Most of the results of this paper were originally pre- 
sented [3(| by one of us (ES) at the Tenth Eastern Grav- 
ity Meeting. Since then, and while completing this paper, 
independent results complementary to those presented 
here have been presented in Ref. [3lj . 
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